FIX: compare to epsilon during rref #90
Conversation
src/structure/matrix.rs
Outdated
| /// | ||
| /// Implementation of [RosettaCode](https://rosettacode.org/wiki/Reduced_row_echelon_form) | ||
| fn rref(&self) -> Matrix { | ||
| let epsilon = 1e-10; |
There was a problem hiding this comment.
This particular value of epsilon feels a bit arbitrary when applied universally here.
Would you be interested in considering floating-point constants like machine epsilon or smallest positive normal number? I also assume the choice of a particular threshold constant should be justified formally if possible
There was a problem hiding this comment.
I agree with @djmaxus, but if we use f64::EPSILON, then since it's too strict, so it may not be practical. How's your opinion? @developing-human @djmaxus.
Furthermore, I think it's better to use relative tolerance as follows:
let max_abs = self.data.iter().fold(0f64, |acc, &x| acc.max(x.abs()));
let epsilon = (max_abs * 1e-12).max(1e-15);There was a problem hiding this comment.
This is good feedback, thank you.
I had been thinking of having the function accept a value for epsilon (sympy does similar) but this would be a breaking change.
I think computing epsilon based on the values in the matrix is a good solution. I confirmed it works for my use case and pushed an update.
|
Happy new year @Axect. Do you need anything else from me to get this merged? |
- Fix numerical instability in RREF (Reduced Row Echelon Form) by comparing to epsilon instead of zero ([#90](#90)) (Thanks to [@developing-human](https://github.com/developing-human)) - Update `pyo3` dependency to 0.27.1 for `plot` feature compatibility ([#89](#89)) (Thanks to [@JSorngard](https://github.com/JSorngard)) - Fix adaptive step size control exponent for embedded Runge-Kutta methods - Add `order()` method to `ButcherTableau` trait for correct exponent `1/(p+1)` - BS23: `1/3`, RKF45/DP45/TSIT45: `1/5`, RKF78: `1/8` - Fix misleading comments on RKF78 BU/BE coefficients
3 participants
I found a matrix that was numerically unstable while calculating reduced row echelon form. I've updated the comparisons to zero so they instead compare to epsilon, which has been more stable for my use case.
I added a test case for the unstable scenario, as well as a more basic test since one wasn't present yet for rref.
Updating numerical crates is outside my usual wheelhouse, so I'm very open to feedback here!